% Code for CC ->PC paired recordings,

% Run this section if you need to average a bunch
fclose('all')
clear all

path1=['C:\Users\Tom�s\Desktop\Temp\'];

fileName=' '

preCh=0; %input presynaptic channel.

dir_to_search=[path1]
formatpattern = fullfile(dir_to_search, ['**/ch' num2str(preCh) '_*.ibw']);
dInfoNS = dir(formatpattern);
[~, idx] =sort([dInfoNS.datenum]); % sort structure indices by date of creation, to avoid 10 coming before 2.

dinfo=dInfoNS(idx); % turn the not sorted structure to a sorted one.

nTrials=numel(dinfo)
%% initialize variables

spikes=[];  %initialize variable for adding presynaptic spike locations in each trial
IPSCAmplitudeBaselineAT=[]; % for all trials
  IPSCAmplitudeStimAT=[]; %for all trials
  isiIPSCBaselineAT=[]; %for all trials
  isiIPSCStimAT=[]; %for all trials

%%



%
%
% fclose('all')
% clear all
%
% path1=['C:\Users\Tom�s\Desktop\Temp\'];
% recNumb=11; %recording number.5




if preCh == 0
    postCh=1;
else
    postCh=0;
end

%% the loop
for trial=1:nTrials
    trial
    dataStruc(1)= IBWread([dinfo(trial).folder '\' dinfo(trial).name]);
    
    postName=dinfo(trial).name;
    postName(3)=num2str(postCh);   

    dataStruc(2)= IBWread([dinfo(trial).folder '\' postName]);
    
    Dx=dataStruc(1).dx;  %take any IBW to find sampling rate.
    Fs=1/Dx;
    
    
    rawTraces(:,1,trial)=dataStruc(1).y;
    rawTraces(:,2,trial)=dataStruc(2).y;
%         rawTraces(:,2,trial)=medfilt1(dataStruc(2).y,20);

    time=([1:numel(rawTraces(:,1))]*Dx)';
    
    %% remove step artifact from presynaptic channel
    
    buffer=100;
    artifact=[1e5 2e5];
    
    rawTraces(:,1,trial)=rawTraces(:,1,trial)-medfilt1(rawTraces(:,1,trial),100); %remove baseline from presynaptic channel, i.e. remove step but keep spikes.
    rawTraces([artifact(1)-buffer: artifact(1) + buffer],1,trial)= NaN;
    rawTraces([artifact(2)-buffer: artifact(2) + buffer],1,trial)= NaN;
    
  %  figure;plot(time,rawTraces(:,1,trial)); %fig to review spike trace
  %  after baseline substraction
    %% plot the cumulative charge in pC
    
    % figure; plot(time,cumsum(rawTraces(:,2)))
    cumCharge=cumtrapz(time,rawTraces(:,2,trial)-mode(rawTraces(:,2,trial)));
    in=figure;
    
    subaxis(2,1,1);
    plot(time,cumCharge,'linewidth',2);
    xlabel('time (s)')
    ylabel('cumulative charge (pC)') % it's cumulative integral
%     title ([num2str(recNumb)])
    title([dinfo(trial).name])
    vline(artifact*Dx,'--r'); %label where the stimulus occurs
    
    
    % Make PSTH
    [pks, locs]=findpeaks(rawTraces(:,1,trial)*-1,'minPeakprominence',30,'minpeakDistance',50); %finds the presynaptic spikes
    subaxis(2,1,2);
    plot(time,rawTraces(:,1,trial))
    
%     
%     A=psth(locs,200,1/Dx,1,numel(rawTraces(:,1,trial)));
%     xlabel('Time (ms)');
%     ylabel('Average firing rate (spikes/s)');
%     
    %% find instantaneous firing rate of presynaptic cell and plot in a separate
    % plot
  
    A=locs(1:end-1);
    B=locs(2:end);
    isi([1:numel(B-A)])=B-A;
    hold on;
    iFiringRate=1./(isi*Dx);
%     figure
%     plot(B*Dx,iFiringRate([1:numel(B)]))
%     xlabel('time (s)')
%     ylabel('Instantaneous firing rate (spikes/s)')
%     vline(artifact*Dx)
% %     
    %% Have a moving window for calculating charge
    
    postBS=[];
    %%postBS=rawTraces(:,2,trial)-mode(rawTraces(:,2,trial)); % this sucks, use 5th percentile instead
    postBS=rawTraces(:,2,trial)-prctile(rawTraces(:,2,trial),5);
%     windMs=1000; %in miliseconds
%     windSam= windMs/(Dx*1e3);
%     for i=1:(numel(postBS)-windSam); % remove window pad on both sides
%         chargeRate(i)=trapz(postBS(i:i+windSam))*Dx;
%     end
%     
%     time2=time(1+windSam/2:end-windSam/2);
%     figure;
%     plot(time2,chargeRate(:));
%     
    
    %% Bin trace for real binned charge without moving window (avoid delay)
    
    binWidth2=0.5; %in seconds
    binWidth2Sam=binWidth2/Dx;
    for i=1:(numel(postBS)-1)/binWidth2Sam
        chargeRate3(i)=trapz(postBS(1+(i-1)*binWidth2Sam : 1+binWidth2Sam+(i-1)*binWidth2Sam))*Dx;
        time3(i)=time(1+(i-1)*binWidth2Sam)+binWidth2/2;
    end
  
    %% find the firing rate for each time bin.
    
    bounds=[5.5 9.5]; %start and end points, 5.5 is 5(start) + windowMs/2 (1s/2))
    timeBin=0.5; % window for calculating firing rate and charge rate.
    spikeTimes=B*Dx; %convert from indices to seconds.
    
    %find the average firing rate and charge rate for each time bin
    nBins=((bounds(2)-bounds(1))/timeBin )+1;
    for i=1:nBins
        midPoint=bounds(1)+(i-1)*timeBin;
        startPoint=midPoint-timeBin;
        endPoint=midPoint+timeBin;
        
        firstSpikeIdx=find(abs(startPoint-spikeTimes)==(min(abs(startPoint-spikeTimes))));
        lastSpikeIdx=find(abs(endPoint-spikeTimes)==(min(abs(endPoint-spikeTimes))));
        meanFR(i)=mean(iFiringRate(firstSpikeIdx:lastSpikeIdx));
        
        % meanChargeRate(i)=mean(chargeRate(((startPoint-0.5)/Dx):((endPoint-0.5)/Dx))); %Substract the timeBin to correct for padded window for calculating chargeRate (i.e. first 500 and last 500 ms are gone)
        
    end
    
   % baselineChargeRate=mean(chargeRate(1:round(4.5/Dx)));
   % percChargeRateReduction(:)=1-(meanChargeRate(:)./baselineChargeRate);
    
%     figure; plot(meanFR(:),percChargeRateReduction(:),'k.','markersize',30)
%     set(gcf,'color','white');
%     set(gca,'FontSize', 18);
%     box off

%% IPSC Analysis section
%find IPSCs (size and location)
[pksIPSC, locsIPSC]=findpeaks(postBS,'minPeakprominence',20,'minpeakDistance',50);

% figure; %figure for validation of IPSC detection
% plot(time,postBS)
% hold on
% plot(locsIPSC*Dx,pksIPSC,'r.','markersize',10)
% for i=1:numel(locsIPSC)
%     plot([locsIPSC(i)*Dx locsIPSC(i)*Dx],[0 pksIPSC(i)]);
% end

%sort into baseline and stim times
locsBaselineIPSC=locsIPSC(locsIPSC<artifact(1)); %find ipscs that ocurr before stim onset
pksBaselineIPSC=pksIPSC(locsIPSC<artifact(1));

locsIPSCStim=locsIPSC(locsIPSC>artifact(1) & locsIPSC<artifact(2)); %ipsc that ocurr between the two artifacts, i.e. during stimulation time
pksIPSCStim=pksIPSC(locsIPSC>artifact(1) & locsIPSC<artifact(2)); %same for peask

%convert IPSC locations into intervals in miliseconds
isiStimIPSC=(locsIPSCStim(2:end)-locsIPSCStim(1:end-1))*Dx*1e3;
isiBaselineIPSC=(locsBaselineIPSC(2:end)-locsBaselineIPSC(1:end-1))*Dx*1e3;


  %% Append data for each loop iteration
  %presynaptic spikes
  spikes=[spikes ; locs]; % stack all spikes together for psth
%  chargeRateAllTrials(:,trial)=chargeRate;
  BinnedFRAllTrials(:,trial)=meanFR;
  %BinnedChargeRateAllTrials(:,trial)=meanChargeRate;
  %fractionalChargeRateReduction(:,trial)=percChargeRateReduction;
%   postSynaptic IPSC
  IPSCAmplitudeBaselineAT=[IPSCAmplitudeBaselineAT ; pksBaselineIPSC];
  IPSCAmplitudeStimAT=[IPSCAmplitudeStimAT ; pksIPSCStim];
  isiIPSCBaselineAT=[isiIPSCBaselineAT ; isiBaselineIPSC]; %names are confusing but i'm just appending each trial to a general variable
  isiIPSCStimAT=[isiIPSCStimAT ; isiStimIPSC]; % same as above.
  realBinnedChargeRate(:,trial)=chargeRate3;
end
    

%% Make PSTH of all trials (presynaptic firing)
A=psth(spikes,200,1/Dx,nTrials,numel(rawTraces(:,1,trial)));
xlabel('Time (ms)');
ylabel('MLI firing rate (spikes/s)');
box off
set(gcf,'color','white');
set(gca,'FontSize', 18);

%plot IPSC charge binned for all trials.
figure;
hold on
for i=1:nTrials
   
plot(time3,realBinnedChargeRate(:,i),'color',[0 0 0 1/4])
end
plot(time3,mean(realBinnedChargeRate,2),'color',rgb('firebrick'),'linewidth',2)
Oylim=ylim;
ylim([0 Oylim(2)])
set(gca,'FontSize', 18);
set(gcf,'color','white');
box off
xlabel('time (s)')
ylabel('IPSC Charge (A.U)')
vline([5 10],'--k');

meanChargeBinned=mean(realBinnedChargeRate,2);
baselineChargeMean=mean(meanChargeBinned(1:10));
ccFiringChargeMean=mean(meanChargeBinned(11:20));
effectSize=(ccFiringChargeMean-baselineChargeMean)/baselineChargeMean 


% 
% 
% % plot IPSC charge rate for all trials (mean)
% figure;plot(time2,mean(chargeRateAllTrials,2),'linewidth',2,'color',[0 0 0 0.9])
% set(gcf,'color','white');
% box off
% set(gca,'FontSize', 18);
% xlabel('time (s)')
% ylabel(' IPSC Charge (pC/s)' )

%plot histograms for IPSC amplitude and isi for stim vs baseline

figure;
histogram(IPSCAmplitudeBaselineAT,'BinWidth', 20)
hold on
histogram(IPSCAmplitudeStimAT,'BinWidth', 20)
set(gca,'FontSize', 18);
set(gcf,'color','white');
box off
ylabel('Events')
xlabel('IPSC amplitude (pA)');
legend


figure;
histogram(isiIPSCBaselineAT,'BinWidth', 5)
hold on
histogram(isiIPSCStimAT,'BinWidth', 5)

set(gca,'FontSize', 18);
set(gcf,'color','white');
box off
xlabel('inter IPSC interval (ms)')
ylabel('Events' )
legend

%plot cumulative distributions for both
% Amplitude
figure
[f,x] = ecdf(IPSCAmplitudeBaselineAT);
if max(x)<700
    f(end+1)=f(end);
    x(end+1)=700;
end
% plot(x,f)
plot(x,f,'color',rgb('deepskyblue'),'linewidth',2)

hold on
f=[];
x=[];
[f,x] = ecdf(IPSCAmplitudeStimAT);
if max(x)<600
    f(end+1)=f(end);
    x(end+1)=600;
end

plot(x,f,'color',rgb('firebrick'),'linewidth',2)
box off
legend('baseline','cc firing')
set(gca,'FontSize', 18);
set(gcf,'color','white');
xlabel('IPSC amplitude (pA)')
ylabel('Cumulative probability')



% ISI distribution
figure
[f,x] = ecdf(isiIPSCBaselineAT);
% if max(x)<700
%     f(end+1)=f(end);
%     x(end+1)=700;
% end
% plot(x,f)
plot(x,f,'color',rgb('deepskyblue'),'linewidth',2)

hold on
f=[];
x=[];
[f,x] = ecdf(isiIPSCStimAT);
% if max(x)<600
%     f(end+1)=f(end);
%     x(end+1)=600;
% end

plot(x,f,'color',rgb('firebrick'),'linewidth',2)
box off
legend('baseline','cc firing')
set(gca,'FontSize', 18);
set(gcf,'color','white');
xlabel('Inter IPSC interval (ms)')
ylabel('Cumulative probability')


% save([path1 '11_21_18_Pair1.mat'], 'spikes','chargeRateAllTrials','BinnedFRAllTrials','BinnedChargeRateAllTrials','fractionalChargeRateReduction' ,'nTrials','dinfo')

